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The previous paper [Veldhorst et al., J. Chem. Phys. 141, 054904 (2014)] demonstrated that the 
isomorph theory explains the scaling properties of a liquid of flexible chains consisting of ten Lennard- 
Jones particles connected by rigid bonds. We here investigate the same model with harmonic bonds. 
The introduction of harmonic bonds almost completely destroys the correlations in the equilibrium 
fluctuations of the potential energy and the virial. According to the isomorph theory, if these 
correlations are strong a system has isomorphs, curves in the phase diagram along which structure, 
dynamics and the excess entropy are invariant. The Lennard-Jones chain liquid with harmonic 
bonds does have curves in the phase diagram along which the structure and dynamics are invariant. 
The excess entropy is not invariant on these curves, which we refer to as “pseudoisomorphs”. In 
particular this means that Rosenfeld’s excess-entropy scaling (the dynamics being a function of 
excess entropy only) does not apply for the Lennard-Jones chain with harmonic bonds. 


I. INTRODUCTION 

The dynamics of a viscous liquid near the glass tran¬ 
sition is very state-point dependent. Relatively small 
changes in temperature or pressure can change the relax¬ 
ation time and the viscosity significantly. The dynamics 
of a liquid are usually dependent on two state variables, 
and this state-point dependence is material specific. A 
general understanding of the state-point dependence of 
the dynamics of liquids has since a long time been a goal 
in field of liquid-state and glass physics [1-3]. 

Experimental results by Tdlle et al. [4, 5] indicated that 
the problem of what controls the relaxation time of vis¬ 
cous liquids might be simplified, if instead of pressure p 
and temperature T, one focuses on density p and temper¬ 
ature. It was found that the dynamics of ort/io-terphenyl, 
measured at different state points can be collapsed onto 
a single curve by plotting it as a function of h{p)/T, with 
h{p) = p'^. Later results have shown that the dynamics of 
many liquids can be collapsed when plotted as a function 
of h[p)/T, albeit with a material-dependent h{p) [6-8]. 
A review by Roland et al. [9] established that for many 
liquids h{p) is well approximated by a power law of the 
density, p '^‘‘, with 7 s being a material specific scaling pa¬ 
rameter. We refer to this as power-law density scaling. 
The fact that the dynamics of some liquids were found 
to be a function of the combined variable p"^"/T indi¬ 
cated that there might be a single underlying quantity 
that “controls” the dynamics. 

Some liquids have strong correlations in the equilib¬ 
rium fluctuations of the energy and pressure [10-12]. 
Specifically, these correlations are found in the config¬ 
urational parts of the energy E and pressure p, i.e., the 
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potential energy U and the virial W, respectively. These 
only depend on the positions of the particles r^, in con¬ 
trast to the kinetic energy K and the temperature, which 
only depend on the momenta of the particles Pi: 

E = K(pi,...,pN)-l-U(ri,...,rjv), (1) 

pV = NkBT(pi,... ,pn) + W(ri,...,rAr). (2) 


The UW correlations are quantified by the standard 
Pearson correlation coefficient 


{AWAU) 
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( 3 ) 


where A denotes the difference from the mean (AC/ = 
U — (U)) and angular brackets denote the NVT en¬ 
semble average (constant number of particles, volume, 
and temperature). Liquids with R ^ 0.9 were initially 
called “strongly correlating”, but since this sometimes led 
to confusion with strongly correlated quantum systems, 
they are now referred to as “Roskilde simple” liquids or 
just “R liquids” [13-21]. 

The discovery of this class of liquids subsequently led to 
the development of the isomorph theory [22, 23]. The iso¬ 
morph theory explains why R liquids have many proper¬ 
ties that make them simpler than other liquids. The main 
prediction is that liquids that belong to this class have 
curves in their phase diagram called isomorphs along 
which many properties are invariant. The invariance in¬ 
cludes the dynamics, structure, and excess entropy (the 
entropy minus the entropy of the ideal gas at same tem¬ 
perature and density) [23]. These liquids thus have a 
phase diagram that is effectively one dimensional for 
many properties (but not for, e.g., the pressure and the 
free energy). 

The isomorph theory provides a theoretical explana¬ 
tion for the empirical power-law density scaling. The 
theory does not assume anything about the functional 
form of h[p), and it was indeed discovered that for many 
model liquids [24-27] and two real liquids [28] that h{p) 
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is not well approximated by a power law if the change 
in density is larger than 10-20%. Power-law density scal¬ 
ing is a good approximation when density changes are 
small, which is often the case in experiments. Simple 
power-law density scaling also works in the case of low- 
density supercritical fluids [29, 30]. Another advantage 
of the isomorph theory is that it provides a prediction 
for the (density-dependent) value of js, which was pre¬ 
viously treated as an empirical parameter. Thus, 7 s can 
be estimated independently of the scaling procedure, al¬ 
though this is much more easily done in computer simu¬ 
lations [23, 27, 31] than in experiments [32, 33]. 

The isomorph theory is consistent with another scal¬ 
ing method, which was proposed by Rosenfeld [34] and 
later in a slightly different form by Dzugutov [35] . In this 
scaling the dynamics of many liquids were found to be a 
function of the excess entropy. Since excess entropy and 
the dynamics are both invariant along isomorphs, this is 
in agreement with the isomorph theory, although the iso¬ 
morph theory does not predict which function expresses 
the relaxation time in terms of the excess entropy. 

Initially, the isomorph theory was tested on sim¬ 
ple atomic model systems [23, 24] and small rigid 
molecules [25]. However, liquids that have been shown 
to obey power-law density scaling are usually organic 
liquids with often internal degrees of freedom. In par¬ 
ticular, many polymeric liquids obey power-law density 
scaling. This led us to investigate the applicability of 
the isomorph theory to flexible chain molecules in our 
previous publication (paper I) [36]. We showed that the 
isomorph theory describes the properties of the flexible 
Lennard-Jones chains very well, except for small devi¬ 
ations due to intramolecular effects. These effects were 
due to the covalent bonds in the chain, which do not scale 
with density. 

In paper I [36], the covalent bonds were simulated 
using a constraint algorithm to keep the bond length 
fixed. This is not the most common way to simulate 
molecules with Molecular Dynamics [37]. In the field of 
(bio)chemistry, for instance, coarse-grained and all-atom 
models generally have force fields that model covalent 
bonds as harmonic springs [38-40] . In this paper we show 
that the bond type has large implications on the applica¬ 
bility of the isomorph theory. We show that the Lennard- 
Jones chains with harmonic bonds have curves in their 
phase diagram along which the dynamics and structure 
are invariant. Using a method that is not dependent 
on the validity of the isomorph theory, we also generate 
curves of constant excess entropy (configurational adia- 
bats). We find that unlike isomorphs, curves of invari¬ 
ant dynamics do not coincide with the configurational 
adiabats. We propose the name “psewdoisomorphs” for 
curves that have invariant dynamics and structure, but 
not constant excess entropy. Our result has important 
consequences for the applicability of Rosenfeld’s excess- 
entropy scaling We find that this scaling does not work 
for Lennard-Jones chain liquids with harmonic bonds. 

The paper is structured as follows. In the next sec¬ 


tion we briefly review relevant aspects of the isomorph 
theory. We then introduce the simulation method and 
the Lennard-Jones chain model in Sec. Ill, where we also 
show how the different bond types affect the dynamics 
and the UW correlations of the liquid. We construct a 
configurational adiabat in Sec. IV and curves of invari¬ 
ant dynamics in Sec. V, and test to which degree these 
curves resemble isomorphs. The findings are summarized 
in Sec. VI. 


II. THE ISOMORPH THEORY 

In the isomorph theory [23] uses reduced units, mak¬ 
ing quantities dimensionless using macroscopic quanti¬ 
ties such as temperature and pressure. Quantities in re¬ 
duced units are denoted by a tilde. Examples are the re¬ 
duced distance f = reduced energy U = U/{kBT), 

and reduced time i = p^^^{kBT with ma be¬ 
ing the average particle mass. Denoting a configuration 
as R = (ri,..., r^r), it is expressed in reduced units as 
(which is equivalent to scaling the configura¬ 
tion to unit density). 

At two state points with densities pi and p 2 pairs of 
configurations exist that have the same coordinates in 
reduced units 


pJ/^Ri =P^/'R2 = R, (4) 


i.e., they are scaled versions of each other. If the two state 
points have temperature Ti and T 2 , they are defined to be 
isomorphic if the Boltzmann factors of all pairs of scaled 
configurations obey [23] 



C 12 exp 


U^\ 
kBT2 ) ’ 


( 5 ) 


with the same constant Ci 2 . In practice, this should 
hold to a good approximation for all physically relevant 
configurations. 

From this definition it follows that the structure of the 
liquid is invariant in reduced units at isomorphic state 
points, because the relative probabilities of configura¬ 
tions that are scaled versions of each other are the same 
at both state points. Another important consequence 
of Eq. (5) is that the excess entropy Sex is the same at 
isomorphic state points [23]. 

Taking the logarithm of Eq. (5) and expressing it in 
reduced units, one finds that 

t/(Ri) = i/(R 2 )-ln(Ci. 2 ). (6) 


In other words, the potential energy landscapes at the 
two state points are simply scaled versions of each other. 
Defining the reduced force as the gradient of the reduced 
potential energy surface F = —VU with V = 
one finds that the forces at two isomorphic state points 
are the same in reduced units, Fi = F 2 , for configu¬ 
rations obeying Eq. (4). A particle’s reduced mass is 
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given by fhi = mi/ma where ma is the average particle 
mass. Using this, Newton’s second law in reduced units, 
rhiii = Fi, leads to invariant dynamics when expressed 
in reduced units [23]. 

Isomorphs can be generated using the property of con¬ 
stant excess entropy. This is done using the scaling ex¬ 
ponent 


(AWAU) 

((AC/)2) 


(7) 


Using the standard fluctuation formulae it 
can be shown that {AWAU) / ((AC/)^) = 

{d (W)/dT)v/{d {U)/dT)v [11]. This can then be 
shown to be equal to the slope of a configurational 
adiabat in the (log T, log p) using the configurational 
version of the V — T Maxwell relation [23] : 


pinr\ 


( 8 ) 


Equations (7) and ( 8 ) can be used to map out a config¬ 
urational adiabat for any system, using he fluctuations 
in U and W. In practice this is done by calculating 7 
from the fluctuations using Eq. (7) and than calculating 
the temperature at another state point with a slightly 
different density using Eq. ( 8 ). 

For a more in depth description of the isomorph the¬ 
ory the reader is referred to a recent feature article [41]. 
It should also be noted that the isomorph theory was 
recently generalized by defining a Roskilde-simple sys¬ 
tem by the property that the order of potential ener¬ 
gies is maintained for uniform scaling of configurations: 
t7(Ra) < U{Rb) ^ U{XRa) < U(ARfc) [42]. For the 
properties considered in this paper the new formulation 
of the theory leads to the same predictions as the original 
formulation of the theory. 


III. MODEL AND SIMULATION METHOD 


A. Simulation method 


We simulated the flexible Lennard-Jones chain (LJC) 
model in the liquid phase. The chains consist of ten par¬ 
ticles. All but the bonded particle pairs in the chains 
interact through the well-known Lennard-Jones potential 


v{r) = 4e 


Y-r''- 

.^x-6- 

Vcj/ 

W/ 


(9) 


cut and shifted at 2.5cr. The interaction between bonded 
particle pairs is modeled by a harmonic spring 


v{r) = 0.5fc(r — cr)^ 


( 10 ) 


with spring constant k = 3000e/(7^. Note the bond 
length is the same as the distance at which the Lennard- 
Jones energy is zero. The simulation time step was 


0.0025, which is sufficiently small for the large spring 
constant used. The simulations were carried out in the 
NVT ensemble, keeping the temperature fixed using a 
Nose-Hoover thermostat. The Nose-Hoover thermostat 
is known to not sample an harmonic potential correctly, 
but this does not affect the results for dense systems like 
the liquids studied here [43]. 

We employed a cubic bounding box with periodic 
boundary conditions containing N = 2000 particles (200 
chains). The simulations were performed using the MD 
code RUMD [44], which is optimized for GPU comput¬ 
ing [45]. In some of the figures below we compare our 
results to simulations of the LJC with rigid bonds. Most 
of the latter data come from paper I, in which the de¬ 
tails of the simulation of the rigid-bond chains can be 
found [36]. 


B. Background of the model 

The Lennard-Jones chain was first simulated by Kre- 
mer and Grest [46-48], who used it as a coarse-grained 
model to study the properties of polymeric liquids. The 
particles in the chain correspond to groups of atoms, like 
one or several CH„ units in an alkane, or one or several 
monomers in a polymer. For this reason the Lennard- 
Jones particles in the chain are referred to as “segments”. 

Starting at the end of the 90’s, extensive simulations 
of the model have been done to investigate the behavior 
of polymer melts around the glass transition [49-54]. At 
that time the model had already undergone some changes 
compared to the original version. The main difference 
was that the new simulations did not cut and shift the 
potential at the minimum, but also included the attrac¬ 
tive part of the Lennard-Jones potential [55], whereas 
the earlier versions cut and shifted the potential at the 
minimum. A second difference was that Kremer and 
Grest [46] did not use Molecular Dynamics, but Langevin 
Dynamics, which includes a stochastic force similar to 
what is done when simulating an implicit solvent [56]. 

In the original version of the LJG model [46], as well 
as in many later simulations where the attractive part of 
the LJ potential is taken into account [49-54] , bonds were 
modeled with the following finitely extensible nonlinear 
elastic (FENE) potential 


^^(0 =-0-5fcr^„^ln 1 


r 


^max 


2 


( 11 ) 


Here, Vmax is the maximum length of the bond at which 
the potential diverges. Since the FENE potential is 
purely attractive, it is used in addition to the Lennard- 
Jones potential. This is in contrast to the harmonic-bond 
chains used in this work, where there is no Lennard-Jones 
interaction between bonded particles. The combination 
of the FENE and the Lennard-Jones potentials results in 
a potential minimum of approximately 0.96 (t. 
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FIG. 1. Comparison of the dynamics of the Lennard-Jones 
chain model with rigid bonds (blue lines) and harmonic bonds 
(red dashed lines). The self-intermediate scattering function 
Fs{q, t) with q = 7.09 of the segments and the center of mass, 
and the autocorrelation function of the end-to-end vector 
(R(O)R(t)) are plotted, (a) At a standard density (p = 1.0) 
the harmonic and rigid bonds give practically the same dy¬ 
namics, while the LJC model with FENE bonds (green dotted 
lines) has faster dynamics due to the smaller bond length, (b) 
At higher density and temperature, the chains with harmonic 
bonds have slightly faster dynamics than the chains with rigid 
bonds. The two state points in (a) and (b) are isomorphic to 
one another for the rigid-bond chains [36]. 

Some more recent simulations of Lennard-Jones chains 
have used harmonic bonds [29, 57-61] and rigid 
bonds [29, 36, 62]. In these studies, the bond length is al¬ 
ways set to a for both the harmonic and the rigid bonds. 
In the case of harmonic springs, the spring constant is 
always k = SOOOe/cr^, which is rather stiff and leads to 
narrow bond length distributions. Therefore the stiff har¬ 
monic bonds and rigid bonds are often considered to be 
equivalent, at least concerning the phase diagram of the 
LJC phase diagram [29, 63]. The phase diagram of the 
LJC model with FENE bonds is not expected to be the 
same due to the different bond lengths. Moreover, the 
shorter bond lengths will have a significant effect on the 
molecular structure. We therefore decided to investigate 
the effect of non-rigid bonds by comparing our previous 
results for rigid bonds [36] with new simulations of the 
LJC model with harmonic bonds. 


C. Effects of bond type 

We compare the dynamics of the LJC models with dif¬ 
ferent bond types in Eig. 1, where we plot the intermedi¬ 
ate scattering function of the segments and the center of 
mass, as well as the autocorrelation function of the end- 
to-end vector (R(O)R(t)). Eigure 1(a) shows the dynam¬ 
ics at a standard dense liquid state point. As mentioned 
earlier [29, 63], the harmonic and rigid bonds give indis¬ 
tinguishable dynamics at this state point. Nevertheless, 
at another state point with higher density and tempera¬ 
ture, the dynamics of the chains with harmonic and rigid 
bonds start to differ, as shown in Eig. 1(b), so the two 
models with different bond types cannot be considered 


FIG. 2. Scatter plots of the potential energy U and virial 
W equilibrium fluctuations for the LJC model with different 
bond types. Where the rigid bond simulations (a) have quite 
strong UW correlations, the use of harmonic bonds (b) de¬ 
stroys the correlations. The FENE bonds (c) have a slightly 
higher correlation coefficient than the harmonic bonds. For 
all three bond types, the data were obtained at the state point 
p= 1.0, T = 0.7. 

to be equivalent for all state points. 

The two state points investigated in Fig. 1(a) and (b) 
have been shown in paper I to be isomorphic to each 
other for the LJC liquid with rigid bonds, i.e., they have 
to a very good approximation the same dynamics and 
(intermolecular) structure. We also showed there that 
the isomorph in the phase diagram was well described by 
the condition 

= Const., (12) 

with h{p) = 2p® °® — determined by empirical scal¬ 
ing [36]. The fact that the LJC model with harmonic 
bonds does not have the same dynamics as the model 
with rigid bonds indicates one of two things: either the 
model with harmonic bonds conforms to the isomorph 
theory but with isomorphs described by a different /i(p), 
or the model with harmonic bonds does not conform to 
the isomorph theory. Earlier it has been shown that the 
LJC liquid with harmonic bonds obeys power-law density 
scaling [29] , which is a good approximation to isomorphic 
scaling in small density ranges. Moreover, the model has 
also been shown to obey Rosenfeld’s excess-entropy scal¬ 
ing [57, 64], which is in agreement with the isomorph 
theory. These facts indicate that the LJC liquid with 
harmonic bonds could be described by the isomorph the¬ 
ory, albeit with an h{p) that is slightly different from the 
rigid bond chains. As we shall see, this is not the case, al¬ 
though the model does have curves of invariant dynamics 
and structure. 

So far it seems like the bond type only has a small 
effect on the behavior of the liquid. This is, however, 
not the case when looking at a prominent property of R 
liquids. Figure 2 shows scatter plots of the fluctuations 
of the potential energy U and the virial W for the three 
different bond types. It is immediately apparent that 
the bond type has a big effect on the correlations of the 
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these two quantities. For the rigid bonds it was already 
shown in paper I that the correlation coefficient of the 
LJC liquid is R = 0.86 at the state point p = 1.0, T = 0.7. 
Fig. 2(b) shows that if the same state point is simulated 
with harmonic bonds, the correlations disappear almost 
completely {R = 0.28), even though the dynamics do not 
change (see Fig. 1(a)). We also show data for the FENE 
bonds in Fig. 2(c). Here the correlations are stronger 
{R. = 0.48) than in the system with harmonic bonds, but 
still far from the value of the chains with rigid bonds. 

Both the strong UW correlations and the existence of 
isomorphs in a liquid’s phase diagram are properties of 
an R liquid, and it has been shown that UW correlations 
imply the existence of isomorphs and vice versa [23]. In 
view of this it is puzzling that earlier investigations of the 
model with harmonic bonds have shown that it obeys 
power-law density scaling [29] and Rosenfeld’s excess- 
entropy scaling [57, 64], indicating that it is a simple 
liquid, while the strong UW correlations are absent, in¬ 
dicating that it is not an R liquid. 


IV. DYNAMICS AND STRUCTURE ALONG A 
CONFIGURATIONAL ADIABAT 

Isomorphs are often created by keeping excess entropy 
constant (Sec. II), and we create a configurational adia- 
bat using the same method. We circumvent the time- 
consuming calculation of the excess entropy by using 
Eq. ( 8 ) to keep Sex constant (without knowing its value). 
We do this by performing an initial equilibrium simula¬ 
tion at the state point {pi,Ti) = (1.0,0.7). 7 is then cal¬ 
culated from the fluctuations in U and W using Eq. (7). 
Note that this is possible even if the liquid does not obey 
the isomorph theory, and the fluctuations in U and W 
are not correlated [11, 23]. We choose a density p 2 for 
the new state point, which is close to the density of the 
initial state point. It is then possible to calculate the 
temperature T 2 at this new density for which the excess 
entropy is identical to the first state point, by rewriting 
Eq. ( 8 ) to T 2 = Ti{p 2 /pi)^■ The procedure is repeated 
several times by doing an equilibrium simulation at the 
new state point to calculate a new value of 7 and find a 
new state point on the adiabat. It is important to choose 
the change in density \p 2 — pi\ small enough because 7 
may change with density. This can be verified by mak¬ 
ing sure that a further decrease of this density difference 
does not have a significant effect on the result. Here, the 
change in density was 0 . 02 , and we obtained a set of state 
points with densities ranging from 0.96 to 1.08. 

Eor R liquids, a configurational adiabat is an isomorph, 
and therefore the dynamics are invariant along a configu¬ 
rational adiabat when plotted in reduced units. We test 
this in Eig. 3 where we plot various dynamical quantities 
of the harmonic bond LJC model for state points on the 
configurational adiabat. The figure contains (a) mean 
square displacements, (b) incoherent intermediate scat¬ 
tering functions Fs(q, t) and autocorrelation functions of 



t 


FIG. 3. Dynamics at four state points with the same excess 
entropy, as a function of reduced time, (a) The mean square 
displacement of the segments and the center of mass of the 
chain, (b) The segmental and center-of-mass incoherent inter¬ 
mediate scattering function Fs{q,t) and the autocorrelation 
of the end-to-end vector (R(O)R(t)). The length of the scat¬ 
tering vector was kept constant in reduced units as q = qp^^^ 
with q = 7.09 approximately at the main peak of the static 
structure factor, (c) The autocorrelation functions of the first 
and the fifth Rouse modes. All these dynamical quantities dif¬ 
fer at the four state points, showing that for the LJC model 
with harmonic bonds, the dynamics are not a function of the 
excess entropy alone. 


the end-to-end vector (R(O)R(t)), and (c) Rouse-mode 
autocorrelation functions (Xp(O)Xp(t)) for p = 1,5. All 
these quantities, which probe the segmental dynamics of 
the individual LJ particles, as well as the chain dynamics, 
are clearly changing along the configurational adiabat. 

The isomorph theory predicts that the structure of the 
liquid is invariant along a configurational adiabat. For 
completeness we also test this with the radial distribu¬ 
tion functions plotted in Fig. 4. The radial distribution 
functions are split in the intermolecular contributions of 
particle pairs in different molecules and intramolecular 
contributions from pairs in the same molecule. The rea¬ 
son is that we have previously shown that only the inter¬ 
molecular contribution is invariant on the isomorph for 
the LJC model with rigid bonds [36]. For the chains with 
harmonic bonds we find qualitatively the same results on 
the configurational adiabat; the average bond length does 
not scale with density, and is therefore not invariant in 
reduced units. The intramolecular structure is thus not 
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FIG. 4. The radial distribution function g{f) along the con¬ 
figurational adiabat, in reduced units. We split the total (seg¬ 
mental) g(f) into (a) intermolecular contributions and (b) in¬ 
tramolecular contributions. For the intermolecular g{f) there 
is a reasonable collapse, but for the intramolecular g{f) there 
are clear differences, especially around the bond length. 


the same at different densities. The intermolecular struc¬ 
ture is reasonably invariant, since the first peak is exactly 
at the same position for the tested densities. However, 
it should be noted the height of the first peak of g{r) 
changes more than what was found for the LJC with 
rigid bonds on the isomorph [36]. 

From the results presented so far we conclude that 
the LJC model with harmonic bonds does not obey the 
isomorph theory or Rosenfeld’s excess-entropy scaling, 
because the dynamics are not invariant on a curve of 
constant excess entropy. This is in disagreement with 
previous results, that have shown that the LJC with 
flexible bonds does obey Rosenfeld’s excess-entropy scal¬ 
ing [57, 64]. The reason for this discrepancy is the way 
that the excess entropy is calculated. The excess entropy 
can been approximated in different ways, the most ex¬ 
act of which is thermodynamic integration. The simplest 
approximation is the pair entropy S 2 which can easily 
be calculated from the radial distribution function, but 
in this case the bonded particle pairs are excluded in the 
calculation of the g{r). A third method that is commonly 
used employs an equation of state developed using Self- 
Associating Fluid Theory [63]. Voyiatzis et al. [64] have 
compared the effect of the different entropy approxima¬ 
tions on the applicability of Rosenfeld’s excess-entropy 
scaling of Lennard-Jones chains with harmonic bonds. 
They found that the scaling works best when S 2 is used, 
ignoring the bonded particle pairs. If the bonded par¬ 
ticle pairs were not removed in the calculation of the 
entropy, as is the case in the thermodynamic integration, 
the transport coefficients could not be collapsed on a sin¬ 
gle curve. Our method of identifying isomorphs using 
Eq. (7) and (8) avoids these problems. 



FIG. 5. Empirical density scaling of the (reduced) relaxation 
time for the harmonic bond chains. On the left the unsealed 
data from the segmental incoherent intermediate scattering 
function Fs(q,i) [q = 7.09p^l/3)) and the autocorrelation 
of the end-to-end vector (R(O)R(t)) are shown for different 
isochores. On the right, the segmental relaxation data have 
been scaled by hand to obtain h{p), and two other measures 
of relaxation times have been scaled by the same factor. 


V. IDENTIFICATION OF A 
PSEUDOISOMORPH 

Galliero et al. [29] have shown that reduced viscosities 
of the LJC model with harmonic bonds can be scaled 
approximately onto a single curve that is a function of 
h(p) = /T. This may seem surprising given our pre¬ 

viously mentioned result that power-law density scaling 
does not work for the LJC model, but in Ref. [29] the 
collapse is not perfect, and the extent of the densities in¬ 
vestigated is not mentioned. Nevertheless, the results of 
Calliero et al. indicate that a curve exists which is simi¬ 
lar to an isomorph, along which dynamics and structure 
are invariant. 

For an R liquid, the curve in the phase diagram de¬ 
scribed by h(p) is called an isomorph. Moreover, not 
only the reduced viscosity, but all dynamical measures, 
the structure in reduced units, and the excess entropy 
are constant on this curve. In this section we construct a 
curve of invariant dynamics and test to what degree it has 
the properties of an isomorph. Since we have shown in 
the previous section that it cannot be a proper isomorph 
because the excess entropy is not constant, we call this 
curve of invariant dynamics a psewdoisomorph. We con¬ 
struct the pseudoisomorph by empirical density scaling 
of the segmental relaxation times, but unlike Ref. [29] do 
not make any assumption about the functional form of 
the scaling function h{p). 

Relaxation times were determined from the incoher¬ 
ent intermediate scattering function of the segments and 
the center of mass, and from the correlation function of 
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FIG. 6 . The values of the harmonic bond /i(p) obtained by 
empirical scaling (blue crosses) compared to the rigid bond 
/i(p) (black dashed line). The data for the harmonic bond 
chains have been fitted to obtain h{p) = (or¬ 

ange dashed line). We included the shape of the configura¬ 
tional adiabat (red crosses) by plotting the temperature at 
each density divided by the temperature at p = 1. The green 
dots show a curve of constant 7 = 7.57. The inset shows 
the same but in a log-log plot, where a constant 7 leads to a 
straight line. 


the end-to-end vector. We defined the relaxation time as 
the time after which the normalized correlation function 
has decayed to 0.2. Unsealed reduced relaxation times of 
the segmental Fg and the end-to-end vector are shown 
in Fig. 5(a) for five isochores with densities ranging from 
p = 0.95 to p = 1.16. It was not possible to go to higher 
relaxation times (lower temperatures) due to crystalliza¬ 
tion at higher densities and phase separation or negative 
pressures at lower densities. Nevertheless a fairly large 
range of temperatures could be reached for some densi¬ 
ties. 

With standard power-law density scaling, a function 
h{p) = p^‘/T is found which collapses the isochoric data 
when plotted versus p^‘, where 7 s is a material-specific 
constant. Instead, we scale each isochore to collapse the 
segmental relaxation times onto the p = 1.00 isochore. 
Thus for each isochore, a scalar h was chosen by hand 
to collapse the segmental relaxation times as functions 
of h/T. The value of the scaling parameter h was found 
independently for each isochore from the segmental relax¬ 
ation times. The results of the scaling in Fig. 5(b) show 
a good collapse for all three measures of the relaxation 
time, even though only the segmental relaxation times 
were used in the scaling procedure. It was found earlier 
with power-law density scaling that both the segmen¬ 
tal and chain dynamics follow the same scaling [65, 66 ], 
and this has also been confirmed for an all-atom polymer 
model [67]. 

The values for h{p) that were obtained from the scaling 


T + harmonic bond 7f from fluctuations 
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FIG. 7. 7 h,(p) for the harmonic bond chains (orange dashed 
line) calculated from the h{p) fit in Fig. 6 . This is compared 
to the rigid bond 'yh{p) (black dashed line) calculated from 
h{p) and the rigid bond 7 values calculated from the UW 
fluctuations (Eq. (7)). For harmonic bonds, the 7 / values 
calculated from the fluctuations (orange crosses) are signifi¬ 
cantly lower. These have been calculated at state points on 
the configurational adiabat. 


are plotted in Fig. 6 (blue crosses). The data for the 
harmonic springs are compared with the fit of h{p) for the 
rigid bond chains (black dashed line) from paper I [36]. 
There is a small but significant difference in the shapes of 
h{p) for the two models. The difference is most obvious 
at high density, where the harmonic bond h{p) is lower 
than the rigid bond h[p). This is in agreement with the 
data in Fig. 1(b), which show that the dynamics of the 
chains with harmonic bonds is faster at high density. To 
keep the dynamics invariant on a pseudoisomorph, the 
temperature on the isomorph of the harmonic bond chain 
should be lower at high densities than for the isomorph 
of the LJC with rigid bonds. Recall that on an isomorph, 
T oc h{p) (see Eq. (12), so also h{p) should be lower for 
the harmonic bond chains. As in the previous paper, we 
have fitted the data to a function of the form h{p) = 
2 / 5 “ — p^. The resulting function h{p) = 2 / 5 ^-^^ — is 
shown as the dashed orange line. 

The inset of Fig. 6 shows the same data in a log-log 
plot. A power-law h{p) = p'^ corresponds to a straight 
line in this plot (green dots). Our data show that this 
is not a good description of the data, which means that 
power-law density scaling is an approximation that only 
works for the smaller density changes (5%), confirming 
previous findings of ours [28]. 

For liquids that obey the isomorph theory, h{p) de¬ 
scribes the shape of configurational adiabats. From 
Eq. ( 8 ) it then follows that 7 as calculated from the fluc¬ 
tuations (from now on denoted by 7 /) is the same as the 
logarithmic derivative of/i(p) [ 68 ]. Eor a pseudoisomorph 
this may not hold, since h{p) does not describe a config¬ 
urational adiabat. It is however still possible to calculate 
the logarithmic derivative of h(p) for the pseudoisomorph 
as 

The result of this is plotted in Eig. 7 (orange dashed line). 














Our values are consistent with the value 7 « 6.5 that 
Galliero et al. found using power-law density scaling. 
Comparing jhip) of the harmonic bonds with jh{p) for 
rigid bonds (dashed black line) we see that 7 has a similar 
magnitude, but a stronger density dependence for the 
chains with harmonic bonds. For rigid bonds ■jhip) and 
7 / are identical. In contrast, the chains with harmonic 
bonds, using Eq. (7) to calculate 7 / from the fluctuations 
(orange crosses), gives 7 / values that are much lower than 
the 7 (p) from the fitted h{p). This confirms that the 
pseudoisomorph is not a configurational adiabat. 

After having established that the pseudoisomorph is 
not a configurational adiabat, we test to which degree 
the former has other isomorph invariants. We obtain 
a set of pseudoisomorphic state points from the fitted 
expression of/i(p) using T = ro(2p"‘'^^ — (Eq. (13)), 
with To = 0.7 the temperature at p = 1. We use densities 
from 0.96 to 1.20, creating a pseudoisomorph along which 
the density changes by 25%. 

In Eig. 8 we plot different dynamical quantities in 
reduced units at the pseudoisomorphic state points. 
These include the segmental and center-of-mass mean 
square displacements and incoherent intermediate scat¬ 
tering functions, and chain specific quantities like the ori¬ 
entational autocorrelation of the end-to-end vector and 
Rouse-mode autocorrelation functions. By definition the 
relaxation times of the segmental intermediate scatter¬ 
ing function are invariant on the pseudoisomorph. The 
data show that the entire shape of the relaxation func¬ 
tions is the same for all quantities. There is only a slight 
deviation at the high Rouse modes, which correspond to 
movements in short subchains. This is identical to what 
was found in paper I for chains with rigid bonds [36] . The 
lower Rouse modes correspond to larger (sub)chains, and 
these are as invariant as the other quantities in the fig¬ 
ure. All correlation functions are found to be much more 
invariant on the pseudoisomorph than on the configura¬ 
tional adiabat (Fig. 3). 

Next we investigate whether the structure of the liq¬ 
uid is also invariant on the pseudoisomorph. We plot 
the radial distribution functions of the pseudoisomorphic 
state points in Fig. 9. As in Fig. 4 the radial distribution 
function is split into an intermolecular contribution (a) 
and an intramolecular contribution (b). As for the rigid- 
bond isomorphs (paper I), we find that the intermolecular 
structure is invariant on the pseudoisomorph, while the 
intramolecular structure is not. The main reason for this 
is the bonded particle pairs in the molecule. 

From Fig. 9(b) it is clear that the behavior of the bonds 
is complicated and varies on the pseudoisomorph, so we 
now investigate this further. For the LJC model with 
rigid bonds, the bonds cannot follow the scaling, because 
their lengths are kept constant in normal units, meaning 
that they change in reduced units. Nevertheless, their be¬ 
havior on the isomorph was rather trivial; they show up 
as delta functions in the g{r) [36]. Fig. 10(a) shows that 
for the harmonic bonds, also the width of the distribution 
changes on the pseudoisomorph. The bond length dis- 



t 


FIG. 8. The dynamics at the pseudoisomorphic state points. 
The pseudoisomorphic state points were chosen to have the 
same relaxation time of the segmental incoherent intermedi¬ 
ate scattering function (procedure described in text), (a) The 
mean square displacement of the segments (solid lines) and 
the centers of mass (dashed lines) (b) The incoherent inter¬ 
mediate scattering function Fa{q,t) of the segments and the 
centers of mass, and the autocorrelation function of the end- 
to-end vector (R(O)R(t)). (c) Some of the autocorrelation 
functions of the Rouse modes (Xp(O)Xp(f)) in reduced units 
for pseudoisomorphic state points. The collapse is good for 
the lower modes, but for the higher modes some deviation 
is seen. The dynamics are invariant on the pseudoisomoph, 
while they ar not at the configurational adiabat (Fig. 3). 


tributions on the isochore and isotherm (Fig. 10(b) and 
(c)) show that the width only depends on temperature, 
as expected from the equipartition theorem. The aver¬ 
age bond length changes slightly on the pseudoisomorph, 
but not enough to be invariant in reduced units. The 
fact that at high temperatures the chains with harmonic 
bonds have faster dynamics than the chains with rigid 
bonds (see Fig. 1(b)) may be related to the wider bond 
length distribution, making it easier for the segments and 
the chain to cross barriers. 

We plot two measures of the molecular size in Fig. 11; 
the mean square end-to-end vector (R^) and the mean 
square radius of gyration {Rg), both in reduced units. 
The size of the chains is not invariant on the pseudoiso¬ 
morph, and the molecular sizes seem to depend only on 
density, since it is almost constant on the isochores. The 
data are very similar to those found for rigid bonds [36] . 
It seems that the effect of temperature is slightly larger 
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FIG. 9. Radial distribution functions g(r) in reduced units on 
the pseudoisomorph. (a) The intermolecular g{r) is invariant 
on the pseudoisomorph. (b) The intramolecular g{r) is not 
invariant. This is mainly due to the bonded particle pairs, 
but there are also significant differences at larger distances. 
The structure is more invariant on the pseudoisomorph than 
on the configurational adiabat (Fig. 4) 



FIG. 10. Bond-length distributions along the pseudoiso¬ 
morph, an isochore, and an isotherm. The width of the 
bond-length distribution changes along the pseudoisomorph 
and isochore, but is constant on an isotherm (equipartition). 
The average bond length changes slightly on the pseudoiso¬ 
morph. Note the that in this figure, the bond length is not 
given in reduced units. 

for the harmonic bonds, which we attribute to the tem¬ 
perature dependence of the the bond length distributions. 


VI. DISCUSSION AND OUTLOOK 

Our analysis of the structure and dynamics shows 
that the LJC model with harmonic bonds has pseudoiso- 
morphs, which are very similar to the isomorphs of the 
LJC model with rigid bonds. Dynamics and structure 
are invariant on these curves, except when very local in¬ 
tramolecular contributions are considered, such as the 



P T 

FIG. 11. The mean square end-to-end vector (R^) and mean 
square radius of gyration (Rg) on the pseudoisomorph, com¬ 
pared to the isotherm and isochore. These measures of the 
molecule size seem only dependent on density, and are not 
invariant on the pseudoisomorph. 


higher Rouse modes and the separation between nearest 
and next-nearest neighbors. The pseudoisomorphs of the 
chains with harmonic bonds have a different shape from 
the isomorphs of the chains with rigid bonds, especially at 
high densities and temperatures. This difference is pre¬ 
sumably caused by the wider bond length distributions 
for flexible bonds at high temperatures. 

The main differences between the pseudoisomorphs in 
this paper and the isomorphs in paper I are related to the 
fluctuations in the energy and the pressure. The UW cor¬ 
relations are weak, and 7 as calculated from the UW fluc¬ 
tuations does not agree with the logarithmic slope of the 
pseudoisomorph in the p, T phase diagram. Therefore the 
excess entropy is not constant on the pseudoisomorph. 
The LJC liquid with harmonic bonds does thus not obey 
the isomorph theory, even though the same model with 
rigid bonds does and has similar dynamics and structure. 

This also means that the liquid with harmonic bonds 
does not obey Rosenfeld’s excess-entropy scaling when 
the bonds are flexible, i.e., the excess entropy does not 
control the relaxation time. This is in disagreement 
earlier results, where the chains with harmonic bonds 
have been shown to obey Rosenfeld’s excess-entropy scal¬ 
ing [57, 64]. The disagreement of these previous results 
with our conclusion may be ascribed to the fact that the 
collapse of the data was not exact in these previous stud¬ 
ies. Moreover, Voyiatzis et al. have shown that excess 
entropy scaling works best when the entropy is approx¬ 
imated by the pair entropy S 2 and the bonded parti¬ 
cle pairs are ignored [64]. We conclude that the flexible 
bonds contribute to the entropy of the system, but this 
contribution is not related to the (long-time) dynamics 
of the system. 

The pseudoisomorphs in this article have been identi- 
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fied using empirical scaling. This is inferior to real iso- 
morphs which can be constructed using Eqs. (7) and (8). 
It would be desirable to be able to find the pseudoiso- 
morphs without reverting to empirical scaling, and work 
is in progress with this aim. 
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